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Abstract 

Background: The second-order, infinite impulse response notch filter is widely used 
to remove electrical power line noise in electrocardiograms (ECGs). However this 
filtering process often introduces spurious ringing artifacts in the vicinity of raw 
signal with sharp transitions. It is challenging to simultaneously remove these two 
types of noise without losing vital information about cardiac activities. 

Objective: Our objective is to devise a method to remove the power-line 
interference without introducing artifacts nor losing vital information. To this end we 
have developed the "hybrid approach" involving two-sided filtration and multi- 
iterative approximation techniques. The two-sided filtration technique can suppress 
the interference but some cardiac components are lost. The lost information can be 
restored using multi-iterative approximation technique. 

Results: For evaluation, four artificial data sets, each including 91 ECGs of different 
heart rates, were generated by a dynamical model. Four publicly-accessible sets of 
clinical data (MIT-BIH Arrhythmia, QT, PTB Diagnostic ECG, and T-Wave Alternans 
Challenge Databases) were also selected. Our new hybrid approach and the existing 
method were tested with these two types of signal under various pre-determined 
conditions. In contrast with the existing method, the hybrid approach can provide 
more than 27.40 dB and 37.77 dB reduction in signal distortion for 95% and 60% of 
artificial ECGs respectively; it can provide in excess of 1 1.78 dB and 17.48 dB 
reduction in distortion for 95% and 60% of these real records respectively. 

Conclusions: Overall, a significant reduction in signal distortion is demonstrated. 
These test results indicate that the newly proposed approach outperforms the 
traditional method assessed on both the artificial and clinical ECGs and suggest it 
could be of practical use for clinicians in the future. 



Background 

Biopotential signals, such as electrocardiogram (ECG), often suffer from power-line 
interference (PLI, 50 or 60 Hz) since the recorded signal is an output of the electric 
fields of coupling states surrounding main power lines (PLs) and the power of the 
body. PLI is probably the most common problem encountered in the processing of 
biopotential signals. Essentially, a notch filter is adapted for minimizing PLI because of 
its ability to reject narrow band noise. Indeed, the second-order, infinite impulse re- 
sponse (IIR) notch filter is routinely applied for this purpose [1-3]. Because of the 
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transient response effect of the notch filter, the impulse response of this type filter gen- 
erally has an oscillatory behavior, which may cause microvolt-level ringing artifacts 
(RAs, typically ranging between 0 and 40 uV) in the immediate regions of input signal 
with sharp transitions. Besides, it will cause undesirable attenuation in signal compo- 
nents at frequencies close to the center frequency (50 or 60 Hz). Tolerable signal dis- 
tortion needs a narrow stopband bandwidth (SBW); however, a narrower SBW results 
in a longer transient response time (TRT); whilst a longer TRT often incurs more ser- 
ious RAs. It is an inherent contradiction. When an ECG signal is being processed, the 
RAs occur in the right side of QRS complexes, and consequently, this implies that 
many cardiac components are lost in ST-T regions. Serious distortion (signal distortion 
caused by the SBW itself and the appreciable RAs) may make the ECG signal more dif- 
ficult to interpret, particularly for the ST-T segment analysis, QT interval estimation, 
the detection of Ventricular Late Potentials (VLPs) and so on [4-6]. Removal of PLI 
however should be done with utmost stringent efforts not to eliminate or distort the 
raw signals without introducing artifacts nor losing vital information [7-9]. Many so- 
phisticated digital methods have been investigated to cope with either 50 or 60 Hz 
interference [10-12], and they satisfy the requirement for suppression and even elimin- 
ation of PLI during ECG signals acquisition. However, it is impossible to design an IIR 
notch filter to remove PLI without causing distortion [13-15], and this problem is still 
unsolved in practice. In this paper we address the challenges to simultaneously remove 
the PLI and RAs without losing critical cardiac components by developing a new 
method which we call the "hybrid approach". 

Being motivated by early pioneering work on investigating the RAs phenomena 
caused by the suppression of PLI, in this paper, the hybrid approach comprising two- 
sided filtration and multi-iterative approximation techniques, is proposed to simultan- 
eously minimizing the PLI and associated RAs. In the first instance, the two-sided fil- 
tration technique is partitioned into four steps, which are applied to eliminate the PLI, 
localize the RAs and remove them afterward, whilst handling the boundary effects 
which are caused by the practical causal filter. To deal with the inherent contradictions, 
next the multi-iterative approximation technique is accomplished in three steps, which 
are adopted to sequentially reconstruct the lost cardiac information. A combination 
may thus prove to be more effective in eliminating these two types of noise. The elab- 
orate scheme of the hybrid approach is stated in Sect. 3. 

From the practical viewpoints of industrial and clinical applications, a bio-model- 
oriented diagnostic signal processing technique should be evaluated on the clinical 
data that are acquired from numerous subjects, as well as the artificial data that 
cover a variety of specific pre-determined conditions. Using artificial ECGs has a typ- 
ical advantage, in that the signal distortion can be precisely calculated, since the ideal 
signal (i.e., "true" signal) can be reset to any desired case. In this study, the perfor- 
mances of the proposed and existing methods are evaluated in detail with artificial 
ECGs which are generated by an open-source program [16,17], as well as four well- 
used, clinical ECG databases. All of these real data are accessible to the public at 
Physionet [18]. We compared these two methods with respect to signal distortion in 
the absence and presence of artificial PLs, respectively. We also examined noise re- 
duction in the presence of artificial PLs. Specifically, these data sets are classified 
into four groups: with and without the addition of artificial electrical PLs under the 
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notch filter with the center frequency of 50 and 60 Hz, respectively. The first two 
groups that without PLs, are designed for quantitatively investigating the signal dis- 
tortion. The other two groups, which mixed with the artificial PLs, are chosen for 
evaluating the distortion in an environment with interference, together with examin- 
ing the capacities of the newly presented and old methods for removing PLI. The re- 
lationship of SBW and TRT of notch filters has not been well delineated, to provide 
more insight in the next section we also detail the key properties of finite impulse re- 
sponse (FIR) and IIR notch filters, so that their properties can be compared with 
each other. Additionally, related RAs are quantitatively examined, and the challenges 
are outlined. 



Problem statement 

A digital notch filter is a band-stop filter that passes all frequency components except 
those lying within a narrow range centered on a center frequency f 0 . The magnitude re- 
sponse of an ideal notch filter may be given as below, 

w)i = {o, nz ( » 

where 0) = 2nflf s is the normalized digital frequency, o) 0 = 2jTf 0 /f s is the normalized cen- 
ter frequency 2itf 0 .f s is the sampling rate, and /is the specified frequency. In practice, 
the notch filter has a SBW at f 0 , that is, Ao = 2jiAf/f s . Ao and Af are normalized and 
digital SBWs, respectively. 

FIR and IIR notch filters 

Let Hj(z) denote the transfer function of a second-order, FIR notch filter, 

H f (z) = l-2yz~ l +z~ 2 (2) 

where y = cos co 0 . Hj(z) is simple and easy to implement. However, a disadvantage in 
using this kind filter is that the SBW of Hy(z) is relatively large, which could not meet 
the specifications [19]. In order to be applicable at narrow SBW situations, a Hj(z) 
based second-order, IIR notch filter is then commonly used [3,20], 

where /? 0 = 1/(1 + A) and A = tan(AW2). We regulate 0 < A < 1 in this study. Eq. (3.a) can 
be rearranged as follows, 

b 0 + b x z~ l + b 2 z~ 2 

HAz) = : - (3.b) 

a 0 + a\Z~ v + a 2 z~ l 

where a 0 = 1, b 0 = b 2 = fio, a 1 = b 1 = - 2yfi 0 , and a 2 = (1 - A)/? 0 . First provided y* - 1, H t 
(z) contains two poles (a x and a 2 ) inside the unit circle \z\ = 1 at z complex plane, 



a ia =p 0 .[y±dr l +\ i -l) (4) 



where a 1 + a 2 = -a lf a x • a 2 = a 2 . The stability and settling time of H^z) are character- 
ized by these two poles. 
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Let x[n] and y[n] be the input and output signals at discrete time n, respectively. This 
filter can be implemented by the following difference equation, 

2 2 

y[n] = a,ky[n-k] + ^ bkx[n-k] (5. a) 

£=1 /c=0 

where the subscript /< refers to the /<th-order index of H^z). By deduction, in essence, 
Eq. (5.a) can be identified as follows, 

y[n\ =J2h[k]x[n-k] (5.b) 
where h[k] represents the impulse response oiH^z) (see Appendix A), 



h[k) 



*o, k = 0 

(-ai)h[k-l] + b u k = l 

(-a^hik-1] + {-a 2 )h[k-2} + b 2 , k = 2 (6) 

(-cn)h[k-l} + (~a 2 )h[k-2], k>3 



Hi(z) is a stable system, since \a 2 \ < 1 and \a^\ < 1 + a 2 . When H^z) contains a pair of 
complex valued poles or a single negative pole, h[k] will cycle back and forth between 
negative and positive during the transient state [21], which indicates that h[k] is associ- 
ated with an oscillatory behavior. 

Recalling Eq. (5.b), let us consider a finite case, this can be written as, 

K 

y[n] = ^h[k]x[n-k] (7) 

where K is a positive integer. One interpretation of Eq. (7) is that it represents a FIR 
notch system. Figure 1(a) -(d) display FIR and IIR notch filters calculated by Eqs. (6) 
and (3.b), respectively. Notably, the higher order K of the filter, the weaker intensity the 
pass -band ripples and greater the selectivity. Consulting Figure 1(a) -(d), it is clear that 
this kind FIR filter has the following limitations: (A) The order K is considerably higher 
than that of an equivalent second-order IIR filter meeting the same requirements. It 
thus has far more computational complexity. (B) Because of many pass-band ripples, 
signals that include the information of interest inside the relevant frequency bands will 
be grossly distorted. This is an issue related to the pseudo-Gibbs phenomena. 



Inherent contradictions 

The Hi(z) is often used for removing PLI. Figure 2 (a) -(c) display an input clinical ECG that 
is corrupted by real PLI, the output of the input signal after passing it to H^z) and the 
relevant residue portion, respectively. In Figure 2(b) we see that PLI is well canceled. 
Ideally, PLI should be eliminated without any undesirable effects to distort the raw signals. 
Unfortunately, this cannot be achieved completely. RAs can sometimes interfere with the 
ST-T regions. The H^z) operates with IIR of oscillation, which probably distorts the input 
as well as attenuates the amplitude by producing RAs. Consequently, important cardiac 
components will be lost. In contrast to the ECG signal, RAs are not remarkable. In other 
words, they may pollute the ST segments (1 to 20 \iV) [8], but they are not distinguishable 
by the naked eye unless the intensity of which is greater than a threshold, such as 1 \iV. 
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Just as in Figure 2(b), it seems as if there were no RAs. However, if we carefully check 
Figure 2(c) (as indicated by arrows), we find that they still exist 

Sharp transitions of signal may generate noticeable and intolerable RAs in the imme- 
diate vicinities of these abrupt changes. Figure 2(d) -(f) is another example, in which we 
clearly see the spurious effect of RAs: RAs with the amplitude up to 30 |iV, as shown in 
Figure 2(f). Therefore, the RAs should be carefully removed to prevent the distortion 
of input signal. The main cause of RAs is due to the abrupt bandstop of Hi(z), spectral 
components that lie within the Af, as well as those close to f 0 ± A/72, will be attenuated; 
this is the frequency-domain description. In the time domain, the cause of RAs is H^z) 
itself: infinite impulse and oscillatory responses. 

In order to quantitatively investigate the relationship between abrupt discontinuities of 
input signal and the corresponded RAs, we begin with the unit impulse signal, which is, 

*[*] = { J' n = Ho (8) 
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Figure 1 Effects of different values of K. (a) FIR filter with K= 50; (b) FIR filter with K= 1 00; (c) FIR filter 
with /C=200; (d) IIR filter with K=°°; (a)-(c) are calculated by Eq. (6), respectively, (d) is calculated by Eq. 
(3.b), which is equivalent to Eq. (5.b). These notch filters are designed at f 5 = 500 Hz with f 0 = 60 Hz, and 
Af = 2.0 Hz. 
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(See figure on previous page.) 

Figure 2 Results of processing of notch filters H,(z). (a) Clinical ECG with 50 Hz real PLI; (b) The output 
of this kind filter when applied to the signal in (a); (c) Differentiated components of (a) and (b); (d) Clinical 
ECG without PLI; (e) The output of the signal in (d) after passing it to this kind filter; (f) Differentiated 
components of (d) and (e); (g) A simulated unit impulse signal (1 mV); (h) The output of H,(z) when applied 
to the signal in (g); (i) Differentiated results of (g) and (h); (j) The distributions of o\ at various SBW; (k) The 
distributions of duration of RAs at various SBW. 



where n 0 is a specific time. For simplifying the mathematics that follows the processing, 
by definition, S[n] starts at 0 (n 0 = 0), and goes to oo. By making S[n] pass through the sys- 
tem Hi(z), and letting yj[n] be the output, v[n] = 8[n] - yfin] the difference. The variance of 
v[n] is then calculated by, 

°l = I>] 2 = E (s[n}-J2h[k}S[n-k]\ 

n =0 n=0 \ k=0 J 

= (l-/z[0]) 2 + ^%] 2 (9.a) 
= 1-2A[0] +J]%] 2 

n=0 



By deduction (see Appendix B), we relate the a 2 and A 



a 2 v =^-X,(^Af) (9.b) 



In most practical applications, provided y 2 + A 2 - 1 < 0 (i.e., tan(AW2) < sin G) 0 , com- 
monly, Aco/2 « G) 0 )> then the pole radius p for H^z) is given by, 



p=\ai\ = \a 2 \ = ^fai 



1+A /.qN 

If y 2 + A 2 -l < 0, Imai = -lma 2 ^ 0 , , A A y J 

v /TIT /n ^"VA/ 

If y = - V 1-A , «! = a 2 < 0 / 

Referring to Eqs. (9.b) and (10), we come to the overall conclusions: (z) If A « 1, a 2 , 
can be loosely interpreted as a linear function of Af. With A/ increasing, more and 
more signal components in the stop and pass bands will be modified in both the ampli- 
tude ("ripple") and the phase, (ii) Eq. (10) expresses that when Af—> 0, then p—>l. In 
other words, the wider the Af, the closer the location of poles to the origin in the z 
plane, meaning the system H t {z) settles more rapidly (also meaning the system has a 
shorter duration of TRT) [22]. Therefore, it indicates that the duration of RAs tends to 
decrease as Af increases. It is an inherent contradiction of notch filters. 

To visualize this problem, next we use Eq. (8) to generate a 10-second-length signal 
that is digitized at 1000 Hz, as shown in Figure 2(g). For convention, let the impulse 
occur at n 0 = 1 s (amplitude of the impulse is 1 mV). By making Af with an increment 
of 0.1 Hz from 0.5 to 10.0 Hz, we calculated the outputs yfin] by Eq. (5.a) at f 0 = 50 Hz 



Zhou and Zhang BioMedical Engineering Online 2013, 12:42 
http://www.biomedical-engineering-online.eom/content/12/1/42 



Page 8 of 24 



with each Af We thus obtained 96 outputs. For each output, at a specified Af the dur- 
ation of RAs is defined as the time from impulse to the point w, at which, 



Figure 2(h) shows the output at Af= 3.0 Hz, and Figure 2(i) displays the correspond- 
ing residual components. Figures 2(j) and (k) plot all the calculated results. Observe 
that these results agree with the aforementioned conclusions. 

The principle of this algorithm 

As we have seen, a cardinal implication is how the QRS complex affects the output of sys- 
tem Hi(z); to the system it always poses a big impulse. Because H^z) is a causal filter, notice- 
ably time-decaying RAs can only occur on the right side of QRS complexes, see Figures 2 
(e), (f ), (h) and (i). This implies that RAs only depend on the input waveforms regardless of 
whether the output waveforms which are of steep transitions or not. In addition, Eqs. (9.b) 
and (10) offer insights on the determinants of H^z), which are uniquely controlled by its Af. 
These key intrinsic properties of H£z) would be applied in this newly developed approach. 

As before, x[n] denotes the input, the number of samples is L and x T [n] denotes the 
counterpart of the signal x[n], that is, x T [n] =x[L -1-n], Then we construct a mirror 
extended signal, 



We explore two-sided filtration and multi- iterative approximation techniques to 
eliminate the probable PLI and RAs which are contained in x me [n\. For the various fil- 
ter parmeters, let y r me [n] denote outputs of the signal x me [n] after passing it to the sys- 
tems Hi(z), y n m ^i\ denote outputs of this new method when applied to x me [n]. 

Two-sided filtration technique 

The design procedure can be summarized as follows: 
Step 1 - Initialization 

Given f andf 0 , at a specified Af, we use Eq. (3.b) to calculate filter coefficients a 2 , b 0 , 

bi and b 2 . bCons = [f/fbl> M represents a round operator,^, = 125 Hz is a constant. 

bCons is an integer that is chosen to accommodate different f. If bCons < 2, letbCons = 2. 

Step 2 - Twice notch filtering 
(z) First notch filter suppressing. We pass the x me [n] to the notch filter, that is the 
Eq. (5. a) with filter coefficients calculated in the previous step, and get the output 
y r me [n] as the system output. Then we can obtain the differential components dy me [n] 
using the derivative filter followed, 



Figure 3(a) illustrates y r me [n] with the lost cardiac components, as indicated by the 
arrows. By means of signal-mirror extension, for a single filtering of x me [n], it 
includes two operations: one is conducted in x[n] in the forward direction, see 
Figure 3(b), which shows the first half of dy me [n] (0 < n < L - 1 ); the other is 




(ii) 




0 < n < L-l 
L < n < 2L-1 



(12) 



d Jme\ n \ = Xme[n]-f me [n] 



(13) 
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(See figure on previous page.) 

Figure 3 Schematic diagram of the two-sided filtration technique, (a) The output y r me [n] of the notch 
filter Hj{z) when applied to the signal in Figure 2(d); (b) First half of dy me [n]; (c) Another half of dy me [n\; (d) 
First half of dy me [n]; (e) Second half of dy me [n]; (f) First half of ls me [n]; (g) Second half of ls me [n]; (h) Segment 
of js me [n\; (i) Residues x me [n]-y^ e [n]; (j) The ultimate output y n me [n] of the signal in Figure 2(d) processed by 
the present method, (c), (e) and (g) show results at relevant positions of the first halves. To this example, it 
is implemented at f 5 = 500 Hz with f Q = 50 Hz, and Af = 3.0 Hz. 



conducted in the counterpart x [n], which is equivalent to be operated in x[n] in the 
backward direction, as displayed in Figure 3(c), which shows the next half of dy me [n] 
(L < n < 2L - 1) at relevant positions of the first half. The "relevant positions" 
means a sample at n is corresponded to the sample at n* (n* = 2L - 1 - n), 
hereinafter the same. Because H t (z) is a causal system, detectable RAs in both halves 
of dy me [n] lie on right and left arms of the same QRS complexes, respectively, 
(zz) Second notch filter suppressing. Sequentially, let dy me [n] be filtered by Eq. (5. a) 
with the same filter, we obtain the output dy me [n]. dy me [n] may involove the PLI, 
RAs and distorted components which are caused by the boundary effect [3]. 
However, dy me [n] can only contains the major information of RAs and distorted 
portions. Figure 3(d) shows the first half of dy me [n] and Figure 3(e) illuminates the 
second half of dy me [n] . In order to localize RAs in the next step, we need the second 
notch filtering operation. In fact, the possible PLI is directly suppressed in this step. 
Step 3 - RAs localization 

Let us first construct a sequence cs me [n] implemented by the derivative operation as below, 



the time delay of Eq. (14) is bCons/2 samples. Next we let cs me [n] pass through a low-pass 
filter shown as follows, and let ls me [n] denote as the relevant output, 

\ _»-lOrder 

H l0 (z) = - (15) 
Y-Z 1 

where lOrder = 4 • bCons, intrinsic delay of Hi 0 (z) is (lOrder - l)/2 samples and the gain is 
lOrder. For extracting the information of RAs, this low-pass filter is introduced to filter out 
fluctuations around the input sample and suppress undesirable outliers. Figure 3(f) displays 
the first half of ls me [n] and Figure 3(g) shows the second half. From Figures 3(f) and (g), it 
is clear that triangle-like spikes localize the RAs. 
Step 4 - Determination of the output 

Likewise, we construct another sequence ds me [n] y which results from the following 
operation, 

ds me [n] = ls me [n]-ls me [n*} (16) 

The use of Eq. (16) makes it possible to distinguish the RAs by means of positive and 
negative samples of ds me [n\. Due to the same consideration in Step 3, we let ds me [n] be 
processed by another low-pass filter defined as below, 




(14) 



i _„-j Order 

= ^pr- 



ay) 
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where j Order = 16 • bCons, the phase delay of Hn(z) is (j Order -1)12 samples and gain is 
jOrder, and let js me [n] denote the output. Figure 3(h) shows the first half of js me [n], and 
we see Figures 3(b) and (c), it is obvious that to each sample at position n: (1) if js me [n] > 
0, it represents that it gathers the information of RAs that lie within one arm of QRS 
complexes (left or right); (2) if js me [n] < 0, it means that it gather the information of RAs 
that lie within the other arm of QRS complexes (right or left), see Figures 3(b), (c) and 
(h). Therefore, to each sample, we eliminate RAs based upon the criteria as below, 

(0 tf>meM < 0, f me [n] = f me [n] + dy me [n\, 

(ii) \ijs me [n] > 0, y n me [n] = f me [rf] + dy me [n*}; 

(Hi) lijs me [n] = 0 and ls me [n] < ls me [n% y n me [n] = f me [n] + dy me [n}; 

(iv) Otherwise, y n me [n] = y r me [rf\ + dy me [n% 

Furthermore, another benefit of criteria (i)-(iv) is avoiding the boundary effect. 

Figure 3(j) shows the output y n m ^\ processed by this new method, and Figure 3(i) 
shows the residue portion x me [n] -y n me [n] . Pay special attention to the details of Figures 3 
(a) and (j), and see the regions with arrows in Figure 3(j), then consult Figures 3(b), (c) 
and (i), we easily find that RAs are mostly eliminated. Note worthily, variables bCons, 
lOrder and jOrder are self-adaptive with respect to f. Although the introduction of 
both-sided filtration increases the memory requirements and the proposed method also 
increases the complexity by introducing logical tests in Step 4, two-sided filtration re- 
stores the phase distortions caused by the one-sided operation. 

Multi-iterative approximation technique 

The single implementation of previous technique has two minor drawbacks, both of 
which are presented in Figure 3(i): (1) it is based upon the assumption that RAs consist 
of non-overlapping portions at positions of both directions of x[n]. In fact, it would not 
be able to distinguish overlapped RAs. However, H t (z) of a small Afmay incur some ap- 
preciable parts (greater than 1 uV) overlapped in the middle of adjacent QRS com- 
plexes, as shown by the second and forth arrows in Figure 3(i). (2) Another major 
limitation encountered in using such a practical filter is that the system will attenuate 
waveforms constituted by high frequencies surrounding f 0 , as we see from positions 
with first and third arrows in Figure 3(i). As in the former disadvantage, intuitively, we 
might apply the H t (z) with a larger Af to overcome this. From Eq (9.b), however, such a 
filter may cause more attenuation in signal components at frequencies close to f 0 , since 
the duration of RAs and the amount of lost components are mutually dependent upon 
each other. 

A novel method that may overcome the fundamental problems is to repeat the pro- 
cessing of two-sided filtration several times with various Af we name it as multi- 
iterative approximation technique. Specifically, it can be partitioned into ternary steps 
as below, 

Step I - Denoising with a large reference Af. In order to make distinctive RAs do not 
overlap each other, the H^z) with a large enough Af should be adapted. Whereas, it 
does not denote that the larger the Af the better the performance. Consulting Eqs. (9.b) 
and (10) and Figures 2(j), (k), it is clear that duration of RAs tends to decrease 
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drastically as Af increases at small Af cases; but it decreases slightly at large Af cases. 
Furthermore, g\ is roughly linear with respect to the Af. Thus, an appropriate 
compromise should be made to fit the task: an empirical and experimental Af= 6.0 Hz 
is employed here. 

Input x me [n] is depicted in Figure 2(d). Figure 4(b) illustrates the residual part x me \vi\-y n me 
[n], which is obtained at Af= 6.0 Hz, and Figure 4(f) shows the related PSD spectrum. In 
the time domain, certain amounts of components of QRS complexes have been obviously 
lost, as we can see from Figure 4(b). It can likewise be seen in Figure 4(f) that parts of the 
signal of interest, which are close tof 0 , have been suppressed because of this large Af in 
the frequency domain. 

Step II - Reconstruction of details with a small target Af In Step I, we achieve a short 
duration of RAs, and thus RAs can be eliminated. However, because of the large SBW 
Af many cardiac components are lost. Thereby, in addition to possible PLs, residual 
part x me [n] -y^e M ^ so contains lost components of original signal. In order to extract 
these useful details, #me M "^me M * s tnen t0 ^ e processed by the two-sided filtration 
technique with a small Af {Af < 6.0 Hz). Let x d me [n] = x me [n] -y n me [n] denote the input, 
y d me [n] the output. 

Step I and Step II are complementary with respect to Eqs (9.b) and (10). To illustrate, 
Figure 4(c) shows the y d me [n] at Af= 2.0 Hz and Figure 4(g) displays the relevant PSD 
spectrum ofy^Jw]. It can be observed from Figures 4(b), (c), (f) and (g) that most 
details have been reconstructed. 

Step III - Reconstruction of slight details with the same target Af We can see the 
vicinities indicated by arrows in Figure 4(g), y d me [n] may still contain a little valuable 
information, just as the results shown in Figure 4(c): the amplitude of which is greater 
than 1 uV. To further reduce the distortion effects, similarly, let x s me [n] = x d me [n] -y d me 
[n], x s me [n] is then be processed by the two-sided filtration technique with the same Af 
in Step II. Let f me [n] as the output. Figures 4(d) and (h) show f me [n] and the relevant 
spectrum, respectively. We see that spectrum content in pass bands of Figure 4(h) is 
flatter than that of Figure 4(g), this represents that most lost details have been 
reconstructed, since the amplitude of f me [n] is substantially less than 1 uV as shown in 
Figure 4(d). 

The ultimate output signal is derived from, 

fnM = x me [n]-K e [n]-f me [n]), (ne [0,1-1]) (18) 

Figure 4(a) displays the residual part x me [n] -y r me [n] obtained by the old method, and 
the relevant PSD spectrum is illustrated in Figure 4(e). Observe that many cardiac 
components, which lie in the pass bands, have been removed, as illustrated at 
positions with the arrow in Figure 4(e). By this new method, lost components have 
been minimized substantially in both time and frequency domains, as shown in 
Figure 4(a) -(d) and Figure 4(e) -(h), respectively. 

Although one should aim to fully restore the lost components, it should be noted at 
this point that the residual coefficients f me [n] may still contain a certain amount of car- 
diac components which fall within this SBW Af (i.e., \f 0 - Af/2,f 0 + A/72]); consequen- 
tially, these components can never be reconstructed due to the so-called "frequency 
overlap", as we can see from Figure 4(d). It is an inherent weakness of H^z). Figure 5 
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Figure 4 Comparison of residual coefficients and relevant PSD spectra, (a) The x me [n 

from the old method (Af=2.0 Hz); (b) Residual part x me [n]-y n me [n] in Step I (Af=6.0 Hz); (c) 



(h) 



200 



250 



-/meM reSU,tS 

Residual part 



y^ e [n] in Step II (Af=2.0 Hz); (d) Residual part y s me [n] in Step /// (Af=2.0 Hz); (e)-(h) show PSD spectra of 
sequences in (a)-(d), respectively. Spectra of (e)-(h) are calculated by Welch's method [26] at f s = 500 Hz. 



shows the flowcharts of the newly proposed algorithm; Figure 5(a) depicts the flowchart 
of two-sided filtration technique; Figure 5(b) represents the flowchart of multi- iterative 
approximation technique. Notably, Af is the only parameter that needs to be specified 
in this method. 



Zhou and Zhang BioMedical Engineering OnLine 2013, 12:42 
http://www.biomedical-engineering-online.eom/content/12/1/42 



Page 14 of 24 



(a) 



Input: x me [n]J s J 0 ,Af 



Step 1: 



Step 2: 



Calculate fli,a2,^0) ^1,^2 using the Eq. (3.^) 



(i) First notch operation: 

VmeW = Fi\ Hi {z) < x me[n] > 

(ii) Second notch operation: 

dy me [n] = x me [n] - y r me [n] ; 

d y'me M = Fi \H,(z) < d yme[n] > 



Step 3: 



Step 4: 







cs me [n] = 
ls me [n] = 1 


dy'me [ n ] ~ dy' me [n — bCons] 

^\h 10 (z) < cs me[n] > 





te [n]=l Sme [n]- 
M = Fi\ Hn (z) 



-ls me [ 
< ds m 



Output yinJyi] is determined by the criteria as below: 



(0 If js me [n] < 0, y n me [n] = y r me [n] + dy me [n] ; 
(ii) If js m e [n] > 0, y n me [n] = y r me [«*] + dy me [n*] ; 
(i/i) If 7'5 me [n] = 0 and fo me [n] < /^[n*], 

(iv) Otherwise, y^n] = y r me [n*] +dy me [n*] 



Output: y^ e [n} 



(b) 



Step I: 



Large reference Af = 6.0 Hz; 
f me [n]=BSFT(x rne [n]J s J 0 ,Af) 



Step II: 



x d me [n}=x me [n]-r me [n] 



With a specified small Af (Af < 6.0 Hz); 
^H=BSFT(4 e [n],/„/ 0) A/) 



Step III: 



x s me [ n ] =x d me [n] -y d me [n\ 



With the same specified small Af; 
^ e H=BSFT« ie H,/„/o,A/) 



The ultimate output: 



)C>]=BSFT(x me [ttU,/o,A/) 

Figure 5 Flowcharts of the newly proposed algorithm, (a) Flowchart of the two-sided filtration 
technique; (b) Flowchart of the multi-iterative approximation technique. Y^ uf [n] = Fi\n( z ) <X in [n]> 
indicates that the input X in [n] is filtered by the system H(z), and the relevant output is F ouf [n]. 
Vmein] = BSFT [x me [n], f s , f 0l Af) indicates that the input x me [n] is processed by the Two-sided filtration 
technique with the specified parameters f sl f 0 and Af, and the relevant output is y%ie[ri\. 



Materials and evaluation 
Artificial and clinical ECG data sets 

We simulated four artificial ECG data sets with f s of 250, 360, 500 and 1000 Hz, re- 
spectively. The generator generates realistic ECGs with user-settable parameters, such 
as "sampling frequency", "internal sampling frequency" and so on. This generator can 
be accessed from [16,17]. For this study, we used the generator to simulate 10-second- 
duration data at "internal sampling frequency" of 2000 Hz (but 720 Hz for ECGs sam- 
pled at 360 Hz) with the specified "mean heart rates", and regulated other parameters 
with default values [17]. To assess the performance of these two methods in an envir- 
onment with various sharp transitions, to each data set, we simulated ECGs with the 
"mean heart rates" ranging from 50 to 140 beats per minute (BPM), in increments of 1 
BPM. We thus obtained 91 data for each data set. 

Four widely used sets of real ECGs (MIT-BIH Arrhythmia Database [MITDB], QT 
Database [QTDB], PTB Diagnostic ECG Database [PTBDB], and the T-Wave Alternans 
Challenge Database [TWADB]) were also selected for evaluation [18], see Table 1 for 
details. All of these data were tested in this study. 
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Performance metrics 

The main power supply is not perfectly stable. In some countries, tolerance of the fre- 
quency variation of PLs is 1% [23]. In poor electrical environments, however, the vari- 
ation may up to 3% [24], Because the bandwidth of real PLs is varying, in order to 
assess the performance of this method under different situations, to both clinical and 
artificial data, we conducted these tests by changing Af of the H^z) from 1.0 to 4.0 Hz 
with an increment of 0.1 Hz, respectively. In addition, to determine the capacity of 
eliminating PLI, artificial PLs that were simulated by sinusoidal functions (with the 
amplitude of 0.1 mV) at industrial frequencies (50 and 60 Hz, respectively), were added 
to the raw data. To each lead, at a specified Af different situations can be categorized 
into four groups: O Atf 0 = 50 Hz, implement H^z) without PLI; © At f 0 = 60 Hz, imple- 
ment Hi(z) without PLI; © At f 0 = 50 Hz, implement H t (z) with the addition of 50 Hz 
PLI; © At Jo = 60 Hz, implement H^z) with the addition of 60 Hz PLI. 

We evaluated the relative distortion produced by two methods. The old method 
means the signals were only filtered by the H^z). To each input signal x[n], we obtained 
the outputs y r [n] and y n [n] by processing input signal x[n] with the IIR filter, that is Eq. 
(5.a) and this new method, respectively. We calculated the relevant ratio of percentage 
root-mean-square difference (rPRD, in units of decibels [dB]), and it is given by, 

a 2 

rPRD = 10logPRD 2 r -10logPRD 2 n = lOlog^ (19) 

G n 

where PRD 2 =^, PRD 2 n =^, Am 2 = %x[n] 2 , a 2 = zj, (*[»]-/ [»]) 2 and a 2 n =f^ 

(x[n]-y n [n]) 2 . rPRD illustrates, for the new method, how distortion of the input x[n] is quan- 
titatively lessened in comparison with the old method. It is not feasible to access the purely 
clinical signals, since these real signals might already contain PLs and other broadband 
noises. We hence let these signals be processed via the new approach first, and let the 
resulting signals be the inputs x[n] here. 

To each data set with a specified group, wherein first, for the results rPRDs calculated 
from all of the various SBWs, we can figure out a threshold, at which, those results 
that, in excess of a certain percentage of total rPRDs, are greater than this threshold. 
Let rPRD \ 95 % and rPRD\ 60 o /o denote the thresholds with percentages of 95% and 60%, 
respectively. To facilitate comparisons with the old method in examining the capability 



Table 1 Four publicly-accessible sets of clinical data are selected for evaluation 



Databases 


«Hz) 


Data numbers/Channels 
(Durations) 


Brief description 


QTDB 


250 


105/2-lead 


It was chosen to represent a wide variety of QRS and ST-T 






(15 min.) 


morphologies. 


MITDB 


360 


48/2-lead 


It was obtained from 47 subjects and contains affluent 






(30 min.) 


arrhythmia information. 


™adb 


500 


100/multi-lead 


Including subjects with risk factors, such as myocardial 






(2 min. (appro.)) 


infarctins, transient ischemia, ventricular and so on, as well as 
healthy subjects. 


PTBDB 


1000 


549/1 5-lead 


It was collected from 290 healthy volunteers and patients. 






(1.92 min.(appro.)) 
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of minimization of distortion at a specific SBW, then, for those results rPRDs outputted 
by the same SBW Af at which (i.e., this single A/), we can define rPRD 

threshold with a percentage of 50%. Therefore, rPRD\ 95 o /o , rPRD\ 60 o /o and rPRD 
resent the overall performances of this newly presented method. 



denote the 



50% re P" 



) corresponding to the histograms in the first four rows, respectively; for a 



Results and discussion 

Test results of the artificial and clinical data sets are shown in Figures 6 and 7, respectively; 
where each histogram illustrates the results of a data set for a specific group. For each artifi- 
cial data set of a specific group, we have 2821 results calculated by Eq. (19). For a specific 
group, we have 6510, 2976, 28892 and 204228 results for QTDB, MITDB, TWADB and 
PTBDB, respectively. Insets, in the last rows of Figures 6 and 7, display the statistical results 

(rPRD 

specific group with varying Af each curve is associated to a histogram in the same 
column. Tables 2 and 3 give the statistical results (rPRD\ 9596 and rPRD\ 6096 ) corresponding 
to the histograms in Figures 6 and 7, respectively. 

Figure 6(a)-(e) and Figure 7(a)-(e) plot the results using the artificial data set and 
QTDB (f s = 250 Hz), respectively. To groups O and ©, as Af increases, it tends to pro- 
duce larger rPRD\^ QO/o for both artificial and clinical data sets as a whole. By contrast, to 

both groups © and ©, as Af increases, it tends to produce rPRD ^ 0/o with relationships 

that look like the exponential decay for artificial data sets, while relationships with 

rPRD^Q 0/o decreasing first and then increasing for QTDB data as a whole. With regard 

to different notch frequencies, the performances are significantly different for both arti- 
ficial and clinical ECG data sets, as we see from Figures 6(e) and 7(e). From Table 2, to 
four groups, the minimum rPRD\ 9596 and rPRD\ 6096 were 28.82 dB and 38.82 dB for 
artificial ECGs, respectively. Similarly, in all of four groups, the minimum rPRD\ 9596 
and rPRD\ 6096 were 15.07 dB and 20.58 dB for QTDB ECGs as shown in Table 3, re- 
spectively. It is worth noting that QTDB data was chosen specifically to contain a broad 
variety of QRS morphologies (i.e., various kinds of abrupt discontinuities) [25]. At this 
point, the results of QTDB ECGs are relatively more objective in the present study. 
Figures 6(f)-(i) and 7(f)-(i) display the results using the artificial data set and MITDB 

(f s = 360 Hz), respectively. Figures 6(j) and 7(j) are the relevant rPRD ^ % relations of 

Figures 6(f)-(i) and 7(f)-(i), respectively. MITDB includes in excess of 53 leads with ab- 
normal, wide size or low intensity QRS waveforms, which are of low frequency compo- 
nents (i.e., slow transitions or smooth variations). Therefore, this turns out that less 
distortion results from the old method, since RAs always occur surrounding the notch 
frequency f 0 , that is the high-frequency end in most ECG spectra. Hence, the perform- 
ance of the new method is not very good for MITDB data, as illustrated in Figure 7(j). 
Likewise, with respect to the different notch frequencies, the performances differ mark- 
edly for both artificial and clinical ECG data sets, as we can see from Figures 6(j) and 7(j). 
However, for the same notch frequency (50 or 60 Hz, respectively), the performances are 
not significantly different for MITDB data without or with the addition of artificial PLs. 
Of these statistical results, as displayed in Tables 2 and 3, for four groups, the minimum 
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SBW(Hz; 

- Without PLI (50 Hz) 



SBW(Hz; 

Without PLI (60 Hz) 



SBW(Hz) 

Additive PLI (50 Hz) 



SBW(Hz) 

- Additive PLI (60 Hz) 



Figure 6 Histograms and relevant statistical results of artificial data sets. From left to right: the four 
columns, (a)-(e), (f)-(j), (k)-(p) and (q)-(u), show results of data sets with f s of 250, 360, 500 and 1000 Hz, 
respectively. For each data set, from top to bottom: the first four rows show results of four groups, 
including groups ©, ©, © and ©, respectively. Four curves in each inset in the last row, which contains 
(e), (j), (p) and (u), correspond to four histograms in the same column and ordered by the legends, 
respectively. 



rPRD\ 9596 and rPRD\ 6096 were 28.91 dB and 38.75 dB for artificial ECGs, respectively; for 
the QTDB recordings with four groups, the minimum rPRD\ 9596 and rPRD\ 6096 were 
11.78 dB and 17.48 dB, respectively. 

The statistical results of artificial data set and TWADB (f s = 500 Hz) are shown in 
Figures 6(k)-(p) and 7(k)-(p), respectively. In terms of artificial ECG data with four 
groups, as Af increases, they exhibit consistent tendencies with these of artificial ECGs 
sampled at 250 and 360 Hz. However, for the TWADB data with four groups, as Af in- 
creases, the performances depict no significant consistency as those of QTDB and 

MITDB, since the rPRD Lq % are divergent at sides of low and large SBWs, while rela- 
tively convergent in the middle of SBWs. In addition, for the QTDB and MITDB data 
sets, results with high probabilities lie in low sides of SBWs (i.e., less than 22 dB), but 
for the TWADB data the results of high probabilities lie within high sides (i.e., greater 
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| 27 \ | 

^''''li''''^' 1234 1234 1234 

SBW(Hz) SBW(Hz) SBW(Hz) SBW(Hz) 



Without PLI (50 Hz) Without PLI (60 Hz) Additive PLI (50 Hz) Additive PLI (60 Hz) 

Figure 7 Histograms and relevant statistical results of real data sets. From left to right: the four 
columns, (a)-(e), (f)-(j), (k)-(p) and (q)-(u), show results of data sets with f s of 250, 360, 500 and 1000 Hz, 
respectively. For each data set, from top to bottom: the first four rows show results of four groups, including 
groups ©, ©, © and ©, respectively. Four curves in each inset in the last row, which contains (e), (j), (p) and 
(u), correspond to four histograms in the same column and ordered by the legends, respectively. 



than 30 dB). In Tables 2 and 3, for the four groups, we see that the minimum rPRD\ 95% 
and rPRD\ 6096 were 27.88 dB and 38.60 dB, 14.66 dB and 24.86 dB, for artificial and 
TWADB data sets, respectively. Furthermore, from Figures 6(e), (j), (p) and Table 3, for 
this clinical data set, the performance of the new method is better than that of the 
QTDB and MITDB data sets as a whole. 



Table 2 Statistical results of all of the four artificial ECG data sets 



f" s (Hz) Without PLI Additive PLI (0.1 mV) 





f 0 = 50 Hz (- 


-O-) 


f 0 = 60 Hz (- 


-©-) 


f 0 = 50 Hz (- 


-©-) 


f 0 = 60 Hz (- 


-©-) 




rPRD\ 95 o /o 
(dB) 


rPRD\ 60 o /o 
(dB) 


rPRD\ 9SO/o 
(dB) 


rPRD\ 60 o /o 
(dB) 


rPRD\ 9SO/o 
(dB) 


rPRD\ 6QO/Q 
(dB) 


rPRD\ 9SO/o 
(dB) 


rPRD\ 6QO/o 
(dB) 


250 


28.82 


38.82 


33.20 


42.53 


29.49 


40.25 


35.93 


45.29 


360 


28.91 


38.75 


34.76 


42.60 


29.67 


40.48 


36.86 


45.60 


500 


28.09 


38.60 


33.01 


41.05 


27.88 


39.20 


34.66 


43.68 


1000 


27.40 


37.77 


32.70 


41.19 


27.62 


38.12 


33.78 


42.69 
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Table 3 Statistical results of all of the four clinical ECG data sets 



f s (Hz) 






Without PLI 






Additive PLI (0.1 mV) 






f 0 = 50 Hz (- 


-O-) 


f 0 = 60 Hz (- 


-©-) 


f 0 = 50 Hz (- 


-©-) 


f 0 = 60 Hz (- 


-O-) 




rPRD\ 9SO/o 
(dB) 


rPRD\ 
(dB) 


60% rP/?D| 95 o /o 
(dB) 


rPRD\ 6QO/o 
(dB) 


rPRD\ 95 o /o 
(dB) 


rPRD\ 60 o /o 
(dB) 


rPRD\ 9SO/o 
(dB) 


rPRD\ 6QO/o 
(dB) 


250 a 


17.05 


21.86 


15.07 


20.58 


17.46 


22.14 


16.02 


21.10 


360 b 


15.25 


19.78 


11.78 


17.48 


15.29 


19.90 


12.24 


17.71 


500 c 


14.66 


26.16 


14.71 


24.86 


15.53 


26.66 


15.76 


26.71 


1000 d 


14.67 


24.20 


15.88 


23.85 


16.58 


25.38 


18.07 


26.07 



J QTDB. b MITDB. C TWADB. d PTBDB. 



It can be clearly seen from Figures 6(q)-(u) and 7(q)-(u) that, for the artificial data set 
and PTBDB (f s = 1000 Hz), respectively, the results of the artificial ECGs reveal no sig- 
nificant difference compared with three earlier used artificial data sets. However, the 

clinical ECGs have significant differences in the histograms as well as the rPRD ^ % re- 
lations to these of three previous clinical ECGs. As far as the histograms of the four 
groups are concerned, they exhibit Gaussian probability distributions, since a huge 
number of results was obtained (a total of 204228 results) for each group, as shown in 

Figure 7(q)-(t). In terms of rPRD ^ 0/o relations, as shown in Figure 7(u), four curves 

show tendencies toward convergence as Af increases. Again from Tables 2 and 3, the 
minimum rPRD\ 9596 and rPRD\ 6096 were 27.40 dB and 37.77 dB for artificial ECGs of 
four groups, respectively; and for the PTBDB recordings, the minimum rPRD\ 9596 and 
rPRD\ 6096 were 14.67 dB and 23.85 dB, respectively. 

Generally, for the artificial ECGs of a specific group, differences of results which were 
obtained from different data sets, show no special significance, since the only differ- 
ences of each data set are the sampling rates f s and the duration L. Regarding the four 
clinical ECG data sets (each including four groups), however, they showed obviously 
different performances, since these data sets were obtained from various subjects and 
with different emphases. In particular, the MITDB recordings, comprising many and a 
broad variety of wide size QRS complexes, show relatively poor performance. From all 
the results of artificial and clinical ECGs sampled at each f s , we can find that the per- 
formance exhibited by the artificial data set is significantly better than that of the rele- 
vant clinical data set. The reason is that the clinical ECGs may contain many low 
frequency leads. In summary, all of these test results indicate that the proposed method 
has the capacity to well reduce the PLI, and simultaneously, greatly minimize the RAs 
for both artificial and clinical ECGs. 



Benefits and limitations 

As previously mentioned, we always want the SBW of a notch filter to be very narrow 
for suppressing PLI. The FIR notch filter, which is defined by Eq. (7), does not encoun- 
ter the "infinite impulse" problem, but it requires a large degree K to meet the specifica- 
tion, and this often comes at the cost of high computational complexity. Additionally, 
another problem with the FIR notch filter is that it may still produce RAs since it is os- 
cillatory in the range of the "finite response", because K is a large number. In contrast, 
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our hybrid approach provides primary benefits in eliminating specific interferences. It 
is not limited to reducing fixed fundamental PLI (50 or 60 Hz) but also applicable to 
removing the high-frequency harmonics for some worse cases, since the f 0 can be tai- 
lored to any desired frequency in this approach. 

Thus far, all of our discussions are based upon the tan(AW2) < sin co 0 . In certain ap- 
plications, however, some rare cases may still exist where tan(AW2) > sin G) 0 , especially 
those ECG monitor systems with low sampling rates f s . In general, such situations 
occur at 7s = 2f 0 + (5y and 0 < « f s . By Eq. (4), consider the following two possible cases, 

(i) y=-l That is,^o = f s /2. According to Eq. (3.a), H^z) represents a first-order, IIR 
notch filter that has only one pole a inside the unit circle, for this case, 

p = \a\ = -a = a 2 = -^r, (a < Q,p<xl/Af) (20) 

1 + A V / 

Likewise, for the input S[n], we can obtain the output variance, 

4=^^(4«*r) (2i) 

Eqs. (20) and (21) demonstrate similar forms with Eqs. (10) and (9.b), respectively; 
it implies that the hybrid approach is also applicable. However, real PLs own a 
frequency bandwidth, the counterparts of PLs that lie within the right side of^o will 
pollute valuable information of low frequencies. It is a limitation not caused by the 
method but f itself. Therefore, we recommend Sf> 4.0 Hz for practical application 
based upon the current study and literature [23,24]. 

(ii) - 1 < y< 0 and y 2 + A 2 - 1 > 0 That is, \a x \ * \a 2 \. This yields (see Appendix C), 



i + Ji- r - 2 + r - 2 A 2 

Max(H, \a 2 \) = \y\. 

(«i < 0,a 2 < 0,Max(|ai|, \a 2 \)*Af) 



(22) 



Max(|tfi|, \a 2 1) denotes the larger one of \ai\ and \a 2 \. For \<xi\ * \a 2 \, the settling 
time of Hi(z) is mainly determined by Maxd^, \a 2 \) [22]. Within this situation, 
using the H^z) of a larger Af would not be able to achieve shorter durations of RAs 
but with more cardiac components lost. The "larger" reference Af is then set to the 
specified target Af in Step I to meet the uniformity within the approach. Thus, it is 
worth emphasizing that more precautions should be taken when Af is adjusted to 
avoiding overlapping RAs in the middle of adjacent QRS complexes. Indeed, it is 
essential to remember this limitation as well, since Eq. (10) is not universal but with 
conditions. 



Conclusions 

Instead of frequency response, specifications of second-order, IIR notch filter may be 
given in terms of the impulse response. In this study, we started with the impulse 
response of H^z), the output <r 2 and the settling time that, concerning its behaviour in 
the time domain, have been quantitatively investigated. Minimizing the RAs involved in 
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the signal of sharp transitions by removing PLI is a great concern in the processing of 
biopotential signal The detection and analysis of the VLPs in the ECG signal is highly 
sensitive to the residual PLI and the RAs after the QRS complexes as a result of the 
filtering technique applied. The artifacts may become considerable in cases of high and 
steep complexes. Owing to the intrinsic properties of Hi(z), we proposed a hybrid ap- 
proach, which utilizes H^z), to reliably suppress PLI as well as to aviod the generation 
of RAs. It is applicable to different f s and easy to implement. In fact, Af is the only par- 
ameter that needs to be specified for this approach. Problems are greatly mitigated via 
these techniques. Sufficient results and performance statistics are provided to validate 
the reliability of this method in the test environment with a variety of conditions (e.g., 
artificial and clinical ECGs, of various f s , with altering Af etc.). An eventual consider- 
ation related to practice is f 0 of PLs. To the artificial PLs used in this study, f 0 was set 
to 50 or 60 Hz without varying. However,^ of real PLs, similar to its bandwidth, may 
fluctuate over a small range. Even so, we can refer to many previous studies for how to 
adaptive tracking of the f 0 with serious drift. 

Appendix A: Proof of the Equation (6) 

A real signal x[n] can be expressed as the sum of linear superposition of unit impulse 
functions S[n - k], 



where, if k = n, 8[n - k] = 1; otherwise, if k * n, 8[n - k] = 0. By definition, to each im- 
pulse function S[n - k], the output is denoted as the impulse response h[n- k\. Consult- 
ing Eq. (5.a), we obtain, 

h[n] = -aih[n-l]-a 2 h[n-2] , . 

+ b 0 S[n] + h6[n-l] + b 2 6[n-2] [ } 

Provide H^z) is a causal system, then we conclude that, 

(i) Jfn = 0,h[n-l]=h[n-2]=0. Ignore S[n - 1] and S[n - 2], since S[n - 1] = S[n - 2] 




(a.l) 



= 0, then, 



h[n] = b 0 



(a.3) 



(ii) If n = 1, because h[n-2] =0, and S[n] = S[n - 2] = 0, which results in, 



h[n] = {-ai)h[n-l] +b x 



(a.4) 



(iii) If n = 2, S[n] = S[n - 1] = 0, we have, 



h\n\ = {-ai)h[n-l] + {-a 2 )h[n-2] + b 2 



(a.5) 



(iv) If n > 3, S[n] = S[n - 1] = S[n - 2] = 0, this yields, 



h[n] = {-a^)h\n-\\ + {-a 2 )h[n-2} 



(a.6) 



We thus obtain Eqs (5.b) and (6) in terms of the filter coefficients a 0 , a lf a 2 , b 0 , b 1 
and b 2 . 
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Appendix B: Proof of the Equation (9.b) 

The v[n] = S[n] -y^n] is equivalent to the input 8[n] to be processed by the following 
comb filter, 

H c {z) = 1-H,(z) = co- - , — 2 (b-1) 

where c 0 =Af> 0 . If - 1, H c (z) also contains two poles at a x and a 2 inside the unit cir- 
cle \z\ at z plane. 
According to Parseval's theorem, o 2 v can be calculated by, 

71 

al = J2v[nr=^-[\H c (^)\ 2 dco 

"=° In (b.2) 

_ 1 / H c (z)H c {z- 1 ) 
2nj J c z 

where f c represents the integral taken around the unit circle in counter-clockwise direc- 
tion, and let 

(b.3) 

Z Z i=l{z-CCi)(\-CCiZ) 



By the Residue theorem, we can obtain Eq. (b.4), 

2 

aJ = ^Res[F(z),a*] (b.4) 

k=l 

By Eq. (b.3), we calculate residue values £, and £ at each pole, 



ffi {a\-a<i)(\-a\a<i) 



(b.5) 



_a2 <X2 {a<i-oi\){\-a\a<i) 
And so, 

Appendix C: Proof related to the Equation (22) 

For convention, we define, 

#-i±5 <«■» 

where, 

? = A/i-r 2 + r 2 A 2 (c.2) 
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Take the derivative of both sides of Eq. (c.2) with respect to A, which results in, 

A<p = y~ 2 (p~ l \k\ (c.3) 
Then, 

(1 +A)Ap-(l + <p)AA 



A0: 



(1 + X)(1 + X + AX) 

_ rVA(i+A)-(i + y ) 

(1+A)(1+A + AA) 1 j 

^(A) 



0(A) 



AA 



where, 



>(A) = r 2 A(i+A)-^(i + ^) 

0(A) = p(l+A)(l+A + AA) 



(c.5) 



First assume the ^"(A) < 0 is true, from Eq. (c.5), this yields, 

r - 2 A(l+A)<^(l + ^) (c.6) 

Solve Ineq. (c.6), we obtain (1 + A) 2 < 0. Due to A > 0, thus the assumption < 0 is 
false, then ^"(A) > 0. Because Q{\) > 0, we therefore have, 



$>***XH V / +A -A/ (c.7) 
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